Import libraries¶
Install them with
pip install -r requirements.txt
In [ ]:
import urbanpy as up
import contextily as cx
import matplotlib.pyplot as plt
import plotly
import plotly.express as px
from mpl_toolkits.axes_grid1 import make_axes_locatable
from tqdm.auto import tqdm
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html from .autonotebook import tqdm as notebook_tqdm
In [ ]:
# Activate tqdm progress bar for pandas.apply
tqdm.pandas()
Get city administrative limits¶
In [ ]:
pos_id = 0 # Result position
riyad = up.download.nominatim_osm("Riyadh Governorate, Saudi Arabia", pos_id) # Nominatim query
In [ ]:
# Plot administrative limits
ax = riyad.plot(facecolor="none", edgecolor="r")
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()
Get city population density¶
In [ ]:
# Query the population data in Saudi Arabia from the repository of Meta Data For Good Population Maps in the Humanitarian Data Exchange platform
population_search_results = up.download.search_hdx_dataset("Saudi Arabia")
In [ ]:
population_search_results
Out[ ]:
| created | name | population | size_mb | url | |
|---|---|---|---|---|---|
| id | |||||
| 1 | 2020-03-05 | sau_general_2020_csv.zip | Overall population density | 45.71 | https://data.humdata.org/dataset/14b288ca-1855... |
| 8 | 2019-09-23 | sau_children_under_five_2020_csv.zip | Children (ages 0-5) | 45.55 | https://data.humdata.org/dataset/14b288ca-1855... |
| 9 | 2019-09-23 | sau_elderly_60_plus_2020_csv.zip | Elderly (ages 60+) | 45.48 | https://data.humdata.org/dataset/14b288ca-1855... |
| 10 | 2019-09-23 | sau_men_2020_csv.zip | Men | 45.66 | https://data.humdata.org/dataset/14b288ca-1855... |
| 11 | 2019-09-23 | sau_women_2020_csv.zip | Women | 45.63 | https://data.humdata.org/dataset/14b288ca-1855... |
| 12 | 2019-09-23 | sau_women_of_reproductive_age_15_49_2020_csv.zip | Women of reproductive age (ages 15-49) | 45.62 | https://data.humdata.org/dataset/14b288ca-1855... |
| 13 | 2019-09-23 | sau_youth_15_24_2020_csv.zip | Youth (ages 15-24) | 45.52 | https://data.humdata.org/dataset/14b288ca-1855... |
In [ ]:
# We selected the index 9 that refers to the Elderly (ages 60+)
saudi_arabia_pop = up.download.get_hdx_dataset(population_search_results, 9)
saudi_arabia_pop.head()
Out[ ]:
| longitude | latitude | sau_elderly_60_plus_2020 | |
|---|---|---|---|
| 0 | 39.847083 | 33.000139 | 0.141545 |
| 1 | 40.316806 | 33.000139 | 0.141545 |
| 2 | 44.159583 | 33.000139 | 0.367002 |
| 3 | 44.166250 | 33.000139 | 0.178008 |
| 4 | 44.197083 | 33.000139 | 0.178008 |
In [ ]:
# Number of population points for the country
saudi_arabia_pop.shape[0]
Out[ ]:
9893484
In [ ]:
# Use the city adm limits to filter the country population data
filtered_pop = up.geom.filter_population(saudi_arabia_pop, riyad)
filtered_pop.head()
Out[ ]:
| longitude | latitude | sau_elderly_60_plus_2020 | geometry | |
|---|---|---|---|---|
| 4699315 | 46.770417 | 25.267083 | 0.035855 | POINT (46.77042 25.26708) |
| 4699316 | 46.805694 | 25.267083 | 0.035855 | POINT (46.80569 25.26708) |
| 4699317 | 46.805972 | 25.267083 | 0.035855 | POINT (46.80597 25.26708) |
| 4699318 | 46.806528 | 25.267083 | 0.035855 | POINT (46.80653 25.26708) |
| 4699319 | 46.806806 | 25.267083 | 0.035855 | POINT (46.80681 25.26708) |
In [ ]:
# Number of population points for the city
filtered_pop.shape[0]
Out[ ]:
529876
In [ ]:
# Plot population points
ax = filtered_pop.plot("sau_elderly_60_plus_2020", markersize=0.01, legend=True)
# Plot administrative limit
riyad.plot(facecolor="none", edgecolor="r", ax=ax)
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead if pd.api.types.is_categorical_dtype(values.dtype):
To improve the interpretability of the population plot we will group them in uniform spatial units known as H3 Hexagons.¶
In [ ]:
# Generate hexagons of resolution 7 (~5.1612km2)
riyad_hexs = up.geom.gen_hexagons(resolution=7, city=riyad)
In [ ]:
# Plot the city H3 hexagons
ax = riyad_hexs.plot(facecolor="none", edgecolor="r")
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()
In [ ]:
# Sum the point's population within each hexagons
merges_hexs = up.geom.merge_shape_hex(
riyad_hexs, filtered_pop, {"sau_elderly_60_plus_2020": "sum"}
) # You can use other aggregation methods like max, min, count, and mean
In [ ]:
# Plot hexagons colored by population density
ax = merges_hexs.plot("sau_elderly_60_plus_2020", legend=True, missing_kwds={"color": "grey"})
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead if pd.api.types.is_categorical_dtype(values.dtype):
In [ ]:
# Get amenities (Point of Interest) from Overpass
gdf_pois, _ = up.download.overpass(type_of_data="node", query={"amenity": ["clinic", "hospital"]}, mask=riyad)
gdf_pois.head()
Out[ ]:
| type | id | lat | lon | tags | geometry | poi_type | |
|---|---|---|---|---|---|---|---|
| 38 | node | 4378219310 | 24.400960 | 46.835173 | {'amenity': 'hospital', 'name': 'مستوصف الحاير... | POINT (46.83517 24.40096) | hospital |
| 7 | node | 1354660020 | 24.583313 | 46.864015 | {'amenity': 'hospital'} | POINT (46.86402 24.58331) | hospital |
| 73 | node | 6953249575 | 24.678829 | 46.776056 | {'addr:city': 'الرياض', 'addr:housenumber': 'ح... | POINT (46.77606 24.67883) | hospital |
| 3 | node | 746341536 | 24.687040 | 46.797248 | {'amenity': 'hospital'} | POINT (46.79725 24.68704) | hospital |
| 74 | node | 7011337045 | 24.694442 | 46.790222 | {'addr:street': 'الزبير بن العوام', 'amenity':... | POINT (46.79022 24.69444) | hospital |
In [ ]:
# Plot points of interests
ax = gdf_pois.plot("poi_type", legend=True, figsize=(10, 10))
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead if pd.api.types.is_categorical_dtype(values.dtype):
In [ ]:
# Start routing server (needs docker)
up.routing.start_osrm_server("gcc-states", "asia", "car") # foot,car,bicycle
Starting server ... osrm_routing_server_asia_gcc-states_car Server was started succesfully
In [ ]:
# Calculate travel times (duration in minutes and distance in km) from hexagons to points of interest
merges_hexs_tt = up.accessibility.travel_times(merges_hexs, gdf_pois)
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/urbanpy/accessibility/accessibility.py:260: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. gdf["lon"] = gdf.geometry.centroid.x /Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/urbanpy/accessibility/accessibility.py:261: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. gdf["lat"] = gdf.geometry.centroid.y 0%| | 0/1627 [00:00<?, ?it/s]
Waiting for server to be ready ...
100%|██████████| 1627/1627 [00:15<00:00, 106.13it/s]
In [ ]:
up.routing.stop_osrm_server("gcc-states", "asia", "car")
Server was stoped succesfully
In [ ]:
merges_hexs_tt.head()
Out[ ]:
| hex | geometry | sau_elderly_60_plus_2020 | lon | lat | nearest_poi_ix | distance_to_nearest_poi | duration_to_nearest_poi | duration_to_nearest_poi_label | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 875355da4ffffff | POLYGON ((46.60462 24.61013, 46.60145 24.59749... | 739.076001 | 46.614256 | 24.601532 | 51 | 8.2799 | 10.111667 | 0-15 |
| 1 | 8753734caffffff | POLYGON ((46.50737 24.92900, 46.50419 24.91635... | 51.435090 | 46.517054 | 24.920377 | 96 | 25.7935 | 22.543333 | 15-30 |
| 2 | 875355891ffffff | POLYGON ((47.00107 24.54897, 46.99787 24.53637... | 1.507293 | 47.010623 | 24.540370 | 1 | 22.5693 | 46.918333 | 45-60 |
| 3 | 875355809ffffff | POLYGON ((47.21802 24.66309, 47.21479 24.65051... | 7.536465 | 47.227535 | 24.654486 | 49 | 65.6733 | 111.828333 | 90-120 |
| 4 | 87535582effffff | POLYGON ((47.17685 24.73132, 47.17362 24.71873... | NaN | 47.186380 | 24.722705 | 49 | 61.4836 | 101.770000 | 90-120 |
In [ ]:
# Plot distance and duration maps
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))
# Plot distance
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=-0.1)
merges_hexs_tt.plot("distance_to_nearest_poi", legend=True, ax=ax1, cax=cax)
cx.add_basemap(ax1, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax1.set_axis_off()
ax1.set_title("Distance to Nearest PoI")
# Plot duration
merges_hexs_tt.plot("duration_to_nearest_poi_label", cmap="magma_r", legend=True, ax=ax2)
# Add a basemap
cx.add_basemap(ax2, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax2.set_axis_off()
ax2.set_title("Duration to Nearest PoI")
plt.show()
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead if pd.api.types.is_categorical_dtype(values.dtype): /Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead if pd.api.types.is_categorical_dtype(values.dtype):
Generate interactive maps¶
In [ ]:
plotly.offline.init_notebook_mode()
In [ ]:
fig = up.plotting.choropleth_map(
merges_hexs_tt, "sau_elderly_60_plus_2020", title="Estimated Elderly (60+) Population - 2020", opacity=0.5
)
# Remove the hexagon outlines to make the map clearer
fig.update_traces(marker_line_width=0)
# Make space for the title
fig.update_layout(margin=dict(l=0, r=0, b=0))
In [ ]:
# Filter out the hexagons without population
merges_hexs_tt_filtered_pop = merges_hexs_tt.query("sau_elderly_60_plus_2020 > 0").reset_index(drop=True)
In [ ]:
# Get ordered category labels
category_orders = merges_hexs_tt_filtered_pop["duration_to_nearest_poi_label"].unique().sort_values()
In [ ]:
fig = up.plotting.choropleth_map(
merges_hexs_tt_filtered_pop,
color_column="duration_to_nearest_poi_label",
color_discrete_sequence=px.colors.sequential.Plasma_r,
category_orders={"duration_to_nearest_poi_label": category_orders},
labels={"duration_to_nearest_poi_label": "Minutes"},
title="Travel Time to Nearest PoI",
opacity=0.5,
)
# Make space for the title
fig.update_layout(margin=dict(l=0, r=0, b=0))
# Remove the hexagon outlines to make the map clearer
fig.update_traces(marker_line_width=0)
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/plotly/express/_core.py:2044: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
In [ ]: